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1  Introduction 

Movies  of  the  Solar  atmosphere  reveal  motions  and  variations  of  brightness. 
In  particular  a  sequence  of  two  coronal  images  may  exhibit  the  plane-of- 
the-sky  component  of  the  speed  combined  with  variation  of  the  signal.  The 
present  work  focuses  on  solar  extreme-ultraviolet  images  as  produced  by  the 
Extreme  ultraviolet  Imaging  Telescope  (EIT)  aboard  the  Solar  Heliospheric 
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Observer  (SoHO).  Our  aim  is  to  estimate  both  the  apparent  motion  vector 
and  the  variation  of  brightness  from  two  successive  images.  The  Motion  and 
brightness  Variation  Tracking  (MoVaTrac)  algorithm  is  based  on  a  multi¬ 
scale  optical  flow  algorithm  derived  from  a  local  gradient-based  technique. 
We  demonstrate  a  new  differential  rotation  measurement  technique  and  the 
identification  of  coronal  events  as  outliers  to  the  differential  rotation  or  as 
regions  exhibiting  a  significant  brightness  variation.  Space  weather  services 
have  motivated  this  study.  The  range  of  potential  interests  includes  but  also 
extends  beyond  early  warnings  of  flares  and  coronal  mass  ejection  (CME)  on¬ 
sets.  It  includes  for  example  studies  of  nanoflares  or  macrospicules,  coronal 
seismology,  MHD  and  EIT  wave  investigations,  etc. 

1.1  New  Challenges  in  Space  Weather  Forecasting 

The  study  of  solar  activity  is  needed  for  the  protection  of  satellites,  human 
space  activities,  etc.  Scientists  investigate  the  influences  of  solar  events  -flares 
and  Coronal  Mass  Ejections  (CMEs)-  on  the  Earth  environment.  One  crucial 
task  in  space  weather  applications  consists  in  predicting  major  geomagnetic 
disturbances  from  the  observation  of  the  Sun  corona  and  in  situ  records  in 
the  interplanetary  space.  Expert  observers  are  carrying  out  the  forecast, 
by  especially  using  extreme  ultraviolet  coronal  images  from  the  Extreme  ul¬ 
traviolet  Imaging  Telescope  (EIT,  Delaboudiniere  et  al.  (1995))  instrument 
aboard  SoHO  (Solar  and  Heliospheric  Observatory).  Early  CME  signatures 
include  intensity  variations  in  EUV  images  such  as  dimmings,  EIT  Waves,  or 
motions  such  as  loop  openings  and  prominence  eruptions.  Since  CMEs  can 
cause  geomagnetic  storms,  a  fast  detection  of  these  signatures  is  very  usefull 
to  the  space  weather  forecasters.  Improving  the  prevision  of  geomagnetic  dis¬ 
turbances  from  solar  observations  is  a  crucial  task  for  space  weather  forecast 
centers.  In  the  low  corona,  CME  onset  signatures  include  filament  eruptions, 
coronal  dimmings,  EIT  waves,  loop  openings,  post-eruption  arcades  observed 
in  the  EUV  and  sigmoid-to-arcade  restructuring  in  soft  X-rays  (Hudson  & 
Cliver,  2001).  Unfortunately,  image  processing  of  solar  EUV  image  sequences 
is  difficult.  The  main  problems  stem  from  the  noise,  the  variations  in  source 
brightness  discussed  in  Sect.  1.2,  the  swift  and  hence  under-sampled  topo¬ 
logical  changes,  the  lack  of  spatial  resolution  (spatial  aliasing)  and  the  trans¬ 
parency. 

1.2  Dynamics  in  the  Corona 

The  dynamic  nature  of  the  solar  atmosphere  has  been  observed  and  modeled 
(Rosner  et  al.,  1978;  Moses  et  al.,  1997;  Schrijver  et  al.,  1999).  The  re- 
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cent  years  have  accentuated  the  irrelevance  of  stasis  in  studying  the  physics 
of  the  chromosphere,  transition  region  and  corona.  Dynamics  is  concerned 
with  the  effects  of  forces  upon  motion.  Reciprocally,  the  knowledge  of  move¬ 
ments  and  accelerations  must  constrain  the  models,  validating  or  falsifying 
consequently  the  underlying  physics.  While  forces  are  invisible,  motion  can 
be  quantitatively  estimated  by  observational  means.  Both  spectroscopy  and 
image  sequences  can  reveal  motion.  Spectroscopy  informs  on  the  line-of-sight 
(LOS)  velocity  component  for  traveling  plasma,  whilst  movies  exhibits  only 
the  transversal  component  of  the  speed.  Both  approaches  have  their  respec¬ 
tive  advantages  and  restrictions,  but  they  are  complementary  when  the  two 
types  of  observations  coincide  (e.g.  Kucera  et  al.,  2003).  The  shifts  must  in 
that  case  be  of  the  order  of  the  resolution  element,  i.e.  macroscopic.  The 
present  work  deals  only  with  the  second  method.  It  specially  focuses  on 
solar  extreme  ultraviolet  (EUV)  images  as  produced  in  surfeit  by  EIT  (De- 
laboudiniere  et  al.,  1995)  on  board  SoHO  (see  for  instance  Fig.  1),  TRACE 
(Handy  et  al.,  1999)  or  SPIRIT-CORONAS  (Oraevsky  &  Sobelman,  2002). 
Yet,  the  presented  developments  has  the  potential  to  be  applied  on  magne¬ 
tographic,  photospheric  and  coronographic  image  sequences  or  future  data 
such  as  SWAP/Proba2,  AIA/SDO  or  SECCHI/STEREO. 

Introducing  the  concept  of  apparent  motion  states  immediately  a  fundamen¬ 
tal  limitation  of  imaging  data:  images  are  projections  in  the  2D  plane  of  the 
sky  of  3D  features.  This  problem  of  perspective  in  the  mostly  transparent 
corona  creates  not  only  geometrical  distortions  of  the  velocity  fields  but  also 
ambiguities  if  more  than  one  “feature”  superimpose  on  the  LOS.  It  also  gen¬ 
erates  annoying  occultations  where  opaque  objects  such  as  prominences  or 
the  solar  disc  itself  are  considered.  However,  apparent  motion  goes  beyond 
mere  projection  effects:  motion  tracking  algorithms,  as  well  as  human  ob¬ 
servers,  can  trace  only  intensity  patterns  but  not  actual  plasma  blobs.  With¬ 
out  spectroscopic  information,  propagating  brightness  patterns  are  a  priori 
indistinguishable  from  actual  motions  (De  Groof  et  al.,  2004).  Fortunately, 
at  the  same  time,  this  ambiguity  gives  access  to  a  range  of  wave  phenomena 
able  to  bring  physical  insights  of  their  own  (e.g.  Deforest  &  Gurman,  1998; 
Klimchuk  et  al.,  2004).  Noglik  et  al.  (2005)  calculate  the  reconnection  rate 
from  the  estimated  velocity  at  which  successive  loops  brighten. 

Nevertheless,  solar  rotation  dominates  the  motion  in  image  sequences.  Re¬ 
cent  works  (e.g.  Brajsa  et  al.,  2001;  Vrsnak  et  al.,  2003)  have  studied  differ¬ 
ential  rotation  in  various  layers  of  the  atmosphere.  The  present  report  uses 
their  results  for  calibration  purposes.  In  a  future  publication,  our  method 
eases  the  survey  of  solar  rotation  over  long  periods.  After  differential  rotation 
compensation,  residual  motion  brings  a  wealth  of  insights.  Coronal  activity 
is  obviously  not  limited  to  mere  topological  evolutions  and  rearrangements. 
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a.  EIT  image  Jj. 


b.  EIT  image  J2- 


Figure  1:  A  sequence  of  two  successive  EIT  images  observed  on  1998-05-03, 
repectively  at  times  21:00:21  (7i)  and  21:12:09  (J2)  in  the  CME  Watch  mode, 
at  the  wavelength  19.5  nm. 

Temperature  and  density  variations  are  ubiquitous.  When  they  combine  with 
instrumental  passbands,  the  observed  volumes  of  plasma  brighten  or  darken. 
As  long  as  too  few  bandpasses  are  recorded,  the  differential  emission  mea¬ 
sure  (DEM)  cannot  be  recovered  and  it  is  necessary  to  separate  motion  from 
brightness  variations  (BV).  We  will  see  in  Sect.  2.3.1  the  required  conditions 
for  such  disentanglement. 

Estimating  both  speed  and  brightness  variations  in  parallel  is  mathematically 
logical  and  also  physically  meaningful;  solar  activity  exhibits  both  motion 
and  brightenings,  flaring  and  dimmings.  A  technique  capable  of  producing 
velocity  and  BV  maps  addresses  essentially  all  phenomena  seen  in  imag¬ 
ing  data.  Space  weather  services  have  motivated  this  study  (e.g.  Hochedez 
et  ah,  2005),  although  the  range  of  potential  interests  extends  beyond  early 
warnings  of  flares  and  Coronal  Mass  Ejection  (CME)  onsets.  It  includes 
for  example  studies  of  nanoflares  or  macrospicules,  coronal  and  seismology, 
MHD  and  EIT  wave  investigations,  etc. 

1.3  EUV  Sequence  Analysis  of  the  Solar  Corona 

There  exist  three  main  classes  of  methods  for  motion  analysis  of  a  sequence  of 
two  images  (Barron  et  ah,  1992)  :  gradient-based  (local  or  global),  matching- 
based  and  feature-based.  Gradient-based  and  matching-based,  such  as  Local 
Correlation  Tracking  (LCT)  are  said  to  be  dense,  contrary  to  feature-based 


methods,  in  the  sense  that  they  provide  an  estimation  of  the  velocity  at  ev¬ 
ery  pixel  and  not  only  at  sparse  features.  Dense  does  not  mean  that  the 
quality  of  the  estimation  is  uniform  over  the  image,  and  we  expect  that  some 
solar  regions  have  better  estimations,  since  motion  is  better  defined  at  these 
locations.  In  feature-based  methods,  the  extraction  of  the  features  is  a  prob¬ 
lematic  preprocessing  step. 

Several  authors  have  proposed  matching  techniques  based  on  the  local  corre¬ 
lation  tracking  (LCT)  for  sequences  of  solar  images  at  the  photospheric  level 
(November  &  Simon,  1988;  Welsch  et  ah,  2004).  This  technique  has  also 
been  modified  by  Roudier  et  al.  (1999)  and  applied  to  TRACE  photospheric 
data  (see  Krijger  et  ah,  2002;  Krijger  &  Roudier,  2003).  In  order  to  process 
vector  magnetograms,  Longcope  (2004)  estimates  the  velocity  by  minimizing 
an  energy  term.  Their  “minimum  energy  fit”  enforces  consistency  with  the 
magnetic  induction  equation.  All  these  techniques  deal  with  photospheric 
velocity  flows. 

To  our  knowledge  and  apart  from  our  early  results  (Gissot  et  ah,  2003),  dense 
velocity  field  estimations,  have  not  been  calculated  from  EUV  sequences  of 
the  corona.  Brajsa  et  ah  (2001)  uses  a  feature-based  method  to  track  bright 
points  and  point-like  structures  (PLS).  Each  feature  is  a  point  in  the  (b,  wro^) 
plane  where  b  is  the  heliospheric  latitude  and  wro ^  is  the  angular  velocity. 
The  authors  then  fit  a  parametric  model  of  differential  rotation  through  the 
cloud  of  feature  points.  This  method  is  sensitive  to  the  reliability  of  the 
feature  extraction.  In  this  report,  we  present  a  multiscale  optical  flow  (OF) 
algorithm  derived  from  the  work  of  Lucas  &:  Kanade  (1981)  (LK),  which 
states  a  local  gradient-based  technique.  Our  aim  is  to  estimate  both  the 
fields  of  apparent  displacement  and  brightness  variation  from  two  successive 
images. 

We  further  demonstrate  a  new  differential  rotation  measurement  and  the 
identification  of  coronal  events  as  outliers  to  the  differential  rotation  or  as 
regions  exhibiting  a  significant  BV.  We  apply  the  algorithm  on  EIT  sequences 
of  May  3,  1998  (Fig.  1)  and  April  17,  1999,  that  have  been  observed  during 
the  CME  Watch  observation  mode  at  the  wavelength  19.5  nm  (Fe  XII  spec¬ 
tral  emission  line).  The  report  is  organized  as  follows  :  in  Sect.  2,  we  present 
the  formulation  of  our  new  method  of  motion  analysis.  In  Sect.  3,  we  show 
the  results  of  the  calibration  of  the  method  on  synthetic  signals  and  discuss 
the  analysis  of  observations  in  Sect.  4. 
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2  MovaTrac:  A  New  Algorithm  Analyzing 
Sequences  of  EUV  coronal  Images 

In  order  to  calculate  the  coordinates  of  the  velocity  vector  for  each  pixel,  we 
can  establish  the  optical  flow  (OF)  using  the  brightness  constancy  assumption 
(BCA).  The  BCA  (Horn  &  Schunck,  1981)  states  that  the  source  of  intensity 
remains  constant  over  time.  This  provides  a  single  equation  per  pixel  -the 
optical  flow  constraint  equation  (OFCE)-  while  there  are  two  unknowns, 
namely  the  2  velocity  vector  coordinates.  This  is  the  aperture  problem.  To 
solve  this  under-determination,  following  Lucas  &  Kanade  (1981),  we  fur¬ 
ther  assume  the  local  uniformity  of  the  velocity  field  in  the  neighborhood  of 
a  pixel.  An  alternative  approach  would  be  to  solve  globally  and  iteratively 
until  convergence  of  the  velocity  field  (Horn  &  Schunck,  1981)  by  adding 
a  constraint  on  its  smoothness.  The  disadvantage  of  this  other  method  is 
that  bad  estimations  can  propagate  to  the  rest  of  the  image  through  the 
smoothing.  Furthermore  a  local  diagnostic  on  the  quality  of  the  estimation 
is  impossible  because  the  solution  is  global.  Other  important  approaches  to 
optical  flow  computation  are  probabilistic  (Simoncelli  et  al.,  1991;  Simoncelli, 
1999)  or  robust  minimization  (Black  &;  Anandan,  1993).  In  our  algorithm,  we 
extend  the  BCA  and  the  OFCE  to  estimate  both  the  velocity  and  the  bright¬ 
ness  variation  (BV)  fields  in  the  special  case  of  coronal  ultraviolet  images.  To 
this  end,  a  multiscale  computation  of  the  flow  iterates  the  estimation  from  an 
upper  scale  to  a  lower  scale,  the  scale  being  the  size  of  the  neighborhood  of 
observation.  The  estimation  is  updated  if  a  quality  index  for  the  estimation, 
that  we  define  in  Sect.  2.3.1,  is  improved  across  scales.  In  this  section,  after 
describing  the  preprocessing  that  we  apply  to  the  images  (Sect.  2.1),  we 
introduce  our  new  symmetric  optical  flow  equations  (SOFA,  see  Sect.  2.2.3). 
It  induces  a  symmetry  constraint  between  the  first  and  the  second  image  and 
we  can  express  the  estimation  in  the  Bayesian  framework. 

2.1  Sequence  Preprocessing 

First  we  remove  the  cosmic  ray  hits  (CRH)  using  a  median  filter  (this  proce¬ 
dure  is  part  of  the  solar  soft  library).  We  apply  a  logarithmic  transform  to 
bring  out  the  low-intensity  part  of  the  signal.  Finally,  we  apply  a  gaussian 
smoothing  to  the  image  before  computation.  This  step  allows  us  to  compute 
robust  spatial  gradient  estimations,  and  is  also  important  to  remove  the  esti¬ 
mation  biases.  Thus,  if  we  note  g  the  gaussian  kernel  of  our  smoothing,  our 
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signal  becomes  Ip  =  ga*  Id-  The  spatial  derivatives  are 


dip 

dx 


8_ 

dx' 

9ga 

dx 

=  <  9i,a>  Id  > 


■(da  *  Id) 

*  Id 


Similarly  we  have  : 


dh 

dy 


=<  92, a,  Id  > 


We  have  access  to  an  analytical  derivative  of  the  preprocessed  signal  I  =  ga*I 
where  a  is  the  scale  of  the  gaussian  function  ga(x)  —  exp  (— %&)•  The 
analytical  expression  of  the  gradient  is  given  by: 


V/  = 


(dl  d[  =  ,d9_ 

dx  ’  dy  dx 


*/, 


2.2  The  optical  flow  equations 

There  are  two  equivalent  formulations  of  the  OFCE:  partial  differential  equa¬ 
tion  (PDE)  and  image  registration.  We  aim  at  estimating  the  deformation, 
if  small  enough,  between  the  first  image  I\  and  the  second  /2.  The  PDE 
formulation  is  detailed  in  appendix  (Sect.  A). 


2.2.1  Image  Registration 

We  note  x  =  (x,y)T  is  the  spatial  position  in  the  image  plane  (plane-of-sky). 
We  can  formulate  the  BCA  as: 

I2(x  +  Sx)  =  h{x).  (1) 

A  linear  approximation  of  this  equation,  using  Taylor  series  at  first  order  on 
the  left  hand  side,  gives: 

£ (x,  Sx)  =  I2(x  +  Sx)  -  1 1  (z) 

-  h(x)  +  V/2  •  Sx  —  Ii(x), 
which  leads  to  a  formula  similar  to  (24): 

V/2  •  Sx  +  I2{x)  -  Ii(x)  —  £(£,  Sx).  (2) 

This  formulation  is  used  by  Lucas  &  Kanade  (1981),  where  the  time  variable 
does  not  appear  explicitely.  The  motion  vector  is  then  defined  as  v  =  Sx/St , 
where  St  is  the  time  distance  between  two  images  7i  and  /2.  This  approach 
leads  to  the  same  equation  as  the  PDE  formulation  (see  Sect.  A),  so  that  in 
practice  they  are  equivalent. 
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2.2.2  Optical  Flow  Computation 

In  the  original  Lucas  &;  Kanade  (1981)  optical  flow  estimation,  in  order  to 
solve  the  aperture  problem,  an  assumption  of  uniformity  is  made  on  the 
velocity  vector:  v  has  to  be  uniform  over  a  local  neighborhood  (see  Fig. 
2).  The  set  of  OFCE  equations  over  this  neighborhood  Q  yields  to  a  linear 
system,  where  the  unknowns  are  the  two  parameters  of  the  deformation  (the 
motion  vector).  The  linear  system  can  be  solved  using  a  weighted  linear 
least-squares  method  (Press  et  al.,  1992). 

The  local  uniformity  assumption  (LUA)  The  local  uniformity  assump¬ 
tion  proposed  by  Lucas  &  Kanade  (1981)  states  that  Eqs.  (2)  or  (24)  should 
be  true  locally  around  the  location  estimation,  that  is  over  a  finite  domain  or 
neighborhood  fl(x,  s )  centered  on  x  and  of  size  (or  scale)  c o.  If  Cl(x)  contains 
N  >  2  pixels,  this  assumption  provides  more  equations  of  type  (2)  than  un¬ 
knowns  8xi  so  that  the  system  can  be  solved. 

The  LK  optical  flow  is  the  least-square  (LS)  minimizer  of  £  in  (2)  on 

SxLS{x)  =  argmin  ||£(5f)|&(f)l  (3) 

Sx 

where  ||£(&c)lln(x)  =  Sili £(£»,&<02  =  £T£  is  the  l 2  vector  norm  of  RN 
corresponding  to  the  finite  neighborhood  Q(f).  This  solution  is  the  maximum 
likelihood  estimator  when  the  residual  term  £  is  a  random  variable  following 
a  gaussian  distribution  7V(0,cr|).  In  practice,  we  assign  coefficients  to  each 
pixel,  or  equivalently  to  each  optical  flow  equation  in  the  neighborhood  fl(x), 
so  that  the  pixel  at  the  center  has  a  greater  statistical  weight  than  at  the 
boundaries  of  Q(x).  For  that  we  use  the  Mahalanobis  distance 

lief  =  evf't, 

which  allows  to  assign  weights  to  each  pixels  xt.  It  is  equivalent  to  minimizing 
||tn  *  £||2  where  *  represents  the  2D  convolution  operator  and  w  is  the 
weighting  function.  In  our  application  w  is  an  isotropic  gaussian  function 
normalized  so  that  ||rc||2  =  1. 

LK  solution  to  optical  flow  We  introduce  the  notation  fx  =  df/dx. 
Solving  Eq.  (3)  is  equivalent  to  finding  the  solution  of  the  linear  system 

ASx  =  b,  (4) 


where: 
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hx.(X})  hy{x  l) 
hx(xN)  hyi^N) 

~(h  - 

-(/a  -  h)(SN) 

The  LS  solution  is: 

SxLS  =  (ArVflA)-lATVflb,  (5) 

where  V$  =  £'[£T£],  is  the  covariance  matrix  of  the  residual  error.  The 
covariance  matrix  of  Sxls  is  given  by: 

Vss=(ATVflA)-1.  (6) 

In  the  rest  of  this  report,  we  will  note  this  covariance  matrix  as  V5~~  and  we 
will  use  it  to  define  our  criterion  for  the  quality  of  the  estimation.  Indeed  the 
diagonal  elements  of  Vg-  are  the  variance  of  the  estimation  parameters.  The 
LS  estimator  converges  to  the  correct  solution  if  the  following  conditions  are 
satisfied:  E[£\  =  0  and  =  cr|/,  plus  an  additional  constraint  on  the  matrix 
A:  det  ATA  ^  0. 

Blank  Wall  and  Aperture  The  blank  wall  effect  occurs  in  ID  when  all 
the  dl/dx  term  vanishes,  i.e.  when  the  signal  gradient  is  null  over  the 
neighborhood.  Then  it  is  not  possible  to  estimate  a  motion  vector,  even 
if  the  object  is  moving.  In  2D,  the  matrix  singularity  det  ATA  =  0  means 
either  that  the  image  has  a  uniformly  vanishing  2D  gradient,  which  is  the  2D 
blank  wall,  or  that  the  image  has  a  unidirectional  anisotropic  structure.  For 
that  reason,  the  matrix  ATA  characterizes  the  gradient  texture  of  the  image 
within  a  given  neighborhood  Q.  Furthermore  one  equation  is  not  enough 
anyway  to  derive  the  two  variables  of  the  motion  vector  5x.  This  is  the 
aperture  problem.  We  need  at  least  one  extra  constraint  to  carry  out  locally 
the  estimation,  namely  in  our  case,  the  local  uniformity  assumption. 

2.2.3  Symmetric  Optical  Flow  Analysis  (SOFA) 

We  opt  for  a  symmetric  formulation  of  the  optical  flow  that  can  involve  a 
prior  information,  or  constraint,  on  the  flow  5x.  To  impose  symmetry  in  the 
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optical  flow  estimation,  we  combine  the  two  reciprocal  equations: 

h{x  +  (<5£  -Sxth))  =  h{x), 

Ii(x  -  (8x  -  8xth))  =  /2(£). 

We  linearize  both  equations  assuming  that  5xth  —  0.  It  is  possible  to  impose 
a  predetermined  8xth ,  for  instance  the  theoretical  differential  rotation.  We 
get: 

V/2-5£+/2(f) -/!(£)  =  6,  (7) 

V/i-tf£+72(£)-/i(£)  =  6-  (8) 

We  note: 

/m  =  II  ^  72  (average), 

Id  —  h  —  h  (difference). 

Adding  and  substracting  (7)  and  (8)),  we  obtain: 

2  x  (v/M  •  8$  +  iD(xj)  =  6  +  6,  (9) 

V/z>-5£  =  el-6.  (10) 

The  goal  of  the  algorithm  is  to  minimize  a  cost  function  ||£||2  now  redefined, 
using  the  parallelogram  identity,  as: 

lien2  =  1(IMI2  +  ini2)  =  ii^~ii2  +  ii^-^ii2.  (ii) 

We  note  8x  =  8xis,Sym  the  vector  that  minimizes  the  quantity  argmin^  ||£(<5£)||n- 
We  can  interpret  this  new  minimization  either  as  a  generalized  Tikhonov  reg¬ 
ularization  or  in  the  Bayesian  framework.  The  advantage  of  this  method  is 
that  the  constraint  (seen  as  a  regularization  constraint  or  a  prior  distribu¬ 
tion)  is  derived  from  the  data  (the  matrix  M)  and  does  not  require  any  extra 
arbitrary  parameter.  Furthermore  the  estimation  is  symmetric:  the  param¬ 
eters  that  minimize  ||£||2  do  not  depend  on  the  ordering  between  I\  and  /2. 

We  note: 

lM,y{x  i) 

I  > 

Im,v{xn)  _ 


A  = 


I  M,x{x  1 ) 
Im,x{xn ) 
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b 


M 


-Jd{x\) 

\  > 
—Id(xn) 

I Dtx{x\)  I D,y{X\ ) 
_  Id,x{Xn)  Id,v{xn ) 


where  Ix  means  dl/dx. 

The  cost  function  ||£||2  can  be  written  as: 


ll«1|2  =  5  (l\AS3-  6||2  +  ||MKf)  =  \  (||AH-  6||2  +  ||«||2mtm)  .  (12) 

In  the  Bayesian  framework,  it  appears  that  imposing  the  symmetry  is  equiva¬ 
lent  to  adding  a  prior  constraint  to  the  vector  8x.  Using  the  Bayes’  theorem, 
we  interpret  Eq.  (11)  as  the  following  conditional  probability: 


Px{8x\ID)  = 


Pt(Id\8x)  x  Pr(&r) 

P^Id) 


(13) 


where 


Pr(/D|&r)  oc  exp 
Pr(<S£)  oc  exp 


(~^(A6x  -  b)TVf;  1(ASx  -  b)^j  , 

^-^5xv(MrVflM)5x^J  . 


(14) 

(15) 


Indeed  Eq.  (10)  can  be  used  to  define  a  prior  distribution  on  5x.  This 
prior  distribution  plays  a  role  when  the  gradient  texture  is  deformed,  and 
the  two  diagonal  elements  MTM  then  represent  the  deformation  rates  along 
the  x  and  y  axes.  When  the  gradient  texture  is  deformed  between  I\  and 
I2,  because  of  occlusion  or  strong  deformation  (due  for  instance  to  heatings), 
then  estimation  is  forced  to  tend  to  the  displacement  vector  used  in  the 
linearization  (here  Xth  is  the  0  vector).  Basically,  it  stabilizes  the  estimation 
towards  the  reference  value  (null  displacement  but  could  be  the  differential 
rotation  for  instance)  when  the  gradient  texture  is  modified  between  7j  and 
/2.  Finally,  we  get  the  solution: 

8x  =  (ATV~lA  +  MTVf1M)~1ArVflb.  (16) 

The  covariance  matrix  of  8x  is  now: 

r&sym  =  +  M'Vf'M)-'.  (17) 
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x  x+2 

estimated  pixels  at  scales  s  and  s+1 
window  Q  (estimation  at  scale  s  and  s+1) 
real  velocity  vector 

outlying  velocity  vector  in  the  LS  estimation  at  scale  s+1 


Figure  2:  The  local  uniformity  assumption  in  the  multiscale  estimation  pro¬ 
cess  with  an  update  between  scale  s  +  1  and  scale  s.  At  the  lower  left  corner, 
the  estimation  at  location  (x  +  2,y  —  2)  is  updated  because  the  quality  Q  has 
been  improved  (see  Sect.  2.3.1). 
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Figure  3:  ID  example  of  an  updating  step  in  the  multiscale  MoVaTrac 
algorithm.  At  scale  s  +  1,  the  estimation  is  updated  from  ( 5xs,8Is )  to 
6IS +i). 

2.2.4  Brightness  Variation  Estimation 

We  extend  the  OFCE  with  51,  that  estimates  the  brightness  variation  be¬ 
tween  the  two  images,  relaxing  the  BCA  assumption.  This  approach  has  been 
suggested  by  Lucas  &  Kanade  (1981)  and  has  been  recently  used  by  several 
authors  (e.g.  Odobez  &  Bouthemy,  1995;  Periaswamy  &  Farid,  2003).  An 
example  of  a  one  dimensional  signal  undergoing  a  translation  plus  a  bright¬ 
ness  variation  is  shown  in  Fig.  3.  The  BV  map  can  be  interpreted  as  the 
variation  of  LOS  emission  parameters  (temperature  and  density). 

The  symmetric  optical  flow  equation  (9)  becomes  (without  the  residual 
term  £): 

V/M  -5x  +  ID +  61  =  0,  (18) 

where  51  —  a  x  cj.  The  quantity  a  is  a  parameter  measuring  the  amount 
of  variation  A I  of  intensity  as  a  function  of  the  pixel  unit  Ax.  It  affects 
the  estimation  covariance  matric  but  not  the  estimation  itself.  It  modifies 
the  dynamics  of  the  estimated  parameter  c/.  From  now  on,  we  will  note 
6  =  (5x,  Sy,  C/)T  the  vector  of  parameters  that  we  want  to  estimate. 
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2.3  The  multiscale  SOFA  computation 

In  order  to  make  the  maps  ( 5x,6y ),  which  corresponds  to  the  velocity  vector 
map,  and  81  (bv)  as  local  as  possible,  we  compute  the  optical  flow  in  a 
multiscale  framework  (see  Fig.  4).  The  scale  is  defined  by  the  scale  parameter 
s  of  the  function  w  (see  Sect.  2.2.2).  The  first  computation  is  carried  out 
at  a  predefined  large  scale  So.  At  this  scale,  a  large  number  of  pixels  are 
used  in  the  estimation.  This  estimation  may  nevertheless  be  strongly  biased. 
Furthermore,  within  a  large  neighborhood  {i.e.  at  large  scale),  the  hypothesis 
made  on  the  optical  flow  (extended  BCA  and  LUA)  are  no  longer  valid.  To 
compensate  for  this  effect,  we  add  a  multiscale  updating  step.  It  consists  in 
an  update  of  the  parameter  estimation,  from  the  upper  scale  s0  to  a  lower 
scale,  if  the  quality  index  Q  of  the  estimation  has  increased.  In  other  terms, 
there  is  an  update  when  the  quality  at  the  smaller  scale  s  +  1  is  higher  than 
the  current  larger  scale  s.  From  our  definition  of  the  quality  index  Q  (see 
Sect.  2.3.1  for  details),  we  know  that  if  the  texture  and  the  similarity  have 
been  improved  from  s  to  s  +  1,  then  there  is  an  update.  We  also  leave  as 
an  option  in  our  implementation  the  propagation  of  the  large  scale  velocity 
to  the  lower  scale,  inducing  some  stiffness  between  features  across  scales, 
i.e.  when  finer  scales  undergo  the  motions  of  large  scales.  Combined  with 
an  multiresolution  pyramid  approach,  it  gives  the  possiblity  to  handle  larger 
motions  (more  than  ~  1.5  pixels). 

2.3.1  Quality  index 

Contrary  to  global  methods  such  as  Horn  k.  Schunck  (1981),  our  optical  flow 
approach  gives  the  possibility  to  carry  out  the  registration  depending  on  the 
neighborhood  content.  For  that  reason,  we  define  an  index  Q  (for  “quality”). 
We  will  use  this  index  for: 

•  interscale  comparison  in  order  to  find  the  best  scale  of  estimation  in  the 
multiscale  implementation  by  comparing  Q  between  successive  scales, 

•  interpixel  comparison  the  pixels  that  have  the  best  estimation  by  or¬ 
dering  the  pixels  according  to  the  Q-criterion. 

In  order  to  illustrate  the  meaning  of  Q,  we  can  use  the  space  of  parameter  9.  If 
we  assume  =  ajl  and  since  ATA  is  symmetric,  we  can  write  ATA  =  PDPT 
with  PT  —  P-1,  which  gives: 

VZ1  =  \PDPT 

Sx 
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Figure  4:  Description  of  the  MoVaTrac  algorithm. 
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We  then  define  our  quality  index  as: 

Q  —  min  eigenvalue  (V^1  j  , 

=  \  x  min(Ai,  A2,A3). 

If  we  note  S  =  and  T  =  min(Ai,  A2,  A3),  we  finally  obtain: 

Q  =  S  x  T.  (19) 


where  S  can  be  seen  as  the  similarity  between  both  images  after  deformation, 
while  T  quantifies  the  amount  of  texture  in  the  observation  window.  Indeed, 


we  have  : 


(20) 


where 

w *  ((J2  +  aSI)(*i  +  Sx)  -  h(Xi))2.  (21) 

*<en 


S~l  is  an  estimation  of  the  quantity  cr|.  As  the  noise  on  images  I\  and  I2 
is  supposed  to  be  gaussian  white  noise  N( 0,  erf),  <x|  can  be  interpreted  as  an 
estimation  of  2a] .  This  estimator  has  well-known  statistical  properties  that 
can  be  used  for  statistical  tests.  It  is  also  bounded  by  the  dynamics  of  the 
signal.  To  check  the  numerical  stability  of  the  matrix  Ar  A,  we  compute  its 
eigenvalues  A  and  look  at  the  minimum  eigenvalue  A min.  In  fact,  this  quantity 
corresponds  to  a  “  texture  ”  criterion.  It  is  homogeneous  with  ]T)  I which 
indicates  that  it  is  high  when  the  texture  is  well  suited  to  the  motion  analysis. 
It  is  low  when  there  is  a  blank  wall  effect  or  a  strong  aperture.  For  this  reason 
we  define  a  “texture”  index  T  equal  to  Xmin  —  min(Ai,  A2,  A3). 


2.3.2  Updating  through  scales 

The  variance  cr|  of  the  residual  is  estimated  by  the  dissimilarity  quantity  a| 
and  the  similarity  is  equal  to  the  inverse  of  S.  Two  conditions  have  to  be 
respected  for  a  reliable  matching  (Shi  &  Tomasi,  1994):  first  the  signal  “tex¬ 
ture”’  (T)  must  be  sufficient  (no  aperture,  that  is  enough  “texture”),  and 
second  the  matching  measured  by  similarity  S  must  be  acceptable  after  the 
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estimated  deformation.  Both  conditions  are  satisfied  when  the  quality  index 
Q  is  high,  where  the  term  “high”  will  be  quantified  in  the  calibration  phase. 
When  T  is  low  (unadapted  texture,  blank  wall),  then  S  can  still  be  high  and 
meaningless.  If  S  is  low,  and  T  is  high,  the  quality  is  low  but  the  estimation 
is  then  indeed  incorrect. 

T  —  A min  is  the  texture  parameter,  and  5  is  a  measure  of  similarity  be¬ 
tween  the  first  image  and  the  second  image  after  warping.  If  T  is  large,  then 
the  quality  is  high  because  the  texture  of  the  signal  enables  motion  analysis 
(no  aperture  problem).  If  the  similarity  is  low,  then  the  quality  is  penal¬ 
ized  because  the  registration  (comparison  of  deformed  image  1^,  warped  into 
reference  frame  of  image  1 1),  and  image  I\  is  bad  (high  dissimilarity  <r|). 

Our  quality  criterion  Q  is  related  to  the  theoretical  standard  deviations 
of  the  parameters  given  by  a  linear  least-squares  fit. 

The  updating  stape  corrects  locally  the  error  on  velocity  estimation  but 
increases  the  variance  of  the  estimation  because  there  are  less  pixels  con¬ 
tributing  to  the  estimation.  The  updating  phase  enables  a  tradeoff  between 
large  scales,  where  the  bias  is  due  to  the  non  uniformity  of  the  deformation, 
and  small-scales,  where  the  variance  of  estimation  is  due  to  the  low  number 
of  pixels  (low  N).  In  our  computations,  the  set  of  scales  is  a  sampling  of  four 
scales.  A  higher  number  of  scales  only  improves  the  estimation  and  removes 
the  threshold  effect  due  to  the  updating  step.  For  a  large  number  of  scales, 
the  parameter  maps  (velocity,  brightness  variations  and  quality  index)  are 
smooth.  The  flowchart  of  the  algorithm  is  shown  in  Fig.  4. 

2.3.3  Error  Analysis 

Our  residual  term  may  not  be  gaussian  for  the  following  reasons: 

•  non-gaussian  noise  implying  a  non  gaussian  residual  distribution  (sta¬ 
tistical  error):  e.g.  poisson  noise,  spatial  noise  :  intrapixel  activity 
caused  by  sub-pixel  motions,  can  be  due  to  transparency  effect  can 
cancel  the  advection  effect  at  small  scales.  For  instance  there  could  be 
subpixel  opposite  motion  of  two  structures,  which  mixture  will  yield  to 
outlying  pixels, 

•  modelling  error  (non-zero  second  or  higher  order  terms):  motions  more 
complex  than  translation,  non  uniform  intensity  variations,  presence  of 
non  negligible  higher  order  terms  in  Eq.  (2)  (systematic  error), 

•  noise  in  the  matrix  A  due  to  noisy  spatial  gradient  estimations. 
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2.4  Postprocessing:  coordinate  transform  rotation  ve¬ 
locity  and  meridional  motions 

After  estimating  the  velocity  parameters  vx  —  5x/5t  and  vy  =  5y/5t,  we 
project  them  on  heliospheric  coordinates,  to  get  the  rotation  velocity  and 
the  meridional  velocity.  In  this  report,  we  only  use  the  synodic  observation. 
The  observed  synodic  rotation  velocity  (observed  from  SoHO  or  from  the 
Earth)  is  corrected  to  get  the  sidereal  velocity,  for  a  comparison  of  the  rota¬ 
tion  velocities  that  we  find.  We  first  transform  the  spatial  coordinates  into 
heliospheric  coordinates. 

3  Calibration 

3.1  Synthetic  images 

We  test  the  efficiency  of  our  algorithm  on  a  couple  of  images.  The  first  image 
is  shown  in  Fig.  5. 

I  -4  -4  |2 

G{x^uauAi)  =  Ai  exp(-1—--^'1-  )  ,  i  =  0,7 

The  parameters  have  been  generated  randomly  in  order  to  produce 

an  image  that  is  similar  to  an  EIT  image. 

In  Fig.  6  we  see  that  above  a  predefined  threshold  q,  the  objects  are  well 
located  and  the  estimated  motion  parameters  are  correctly  estimated  (same 
for  bv,  not  show  here). 

3.2  Semi-artificial  EIT  sequences 

We  calibrate  our  method  on  semi-artificial  solar  sequences.  We  apply  a  dis¬ 
placement  field  to  an  EIT  image  I\  and  we  obtain  a  second  image  / i)(*r  using 
a  parametric  model  of  differential  rotation  estimated  by  Vrsnak  et  al.  (2003). 
The  interpolation  procedure  is  based  on  the  IDL  function  of  cubic  interpola¬ 
tion  (Park  &  Schowengerdt,  1983).  This  method  uses  cubic  polynomials  to 
approximate  the  optimal  sine  interpolation  function.  Due  to  spatial  aliasing, 
this  method  gives  bad  results  when  there  are  high  spatial  frequencies  since  it 
assumes  that  the  signal  is  band-limited  (DeForcst,  2004).  This  effect  lowers 
the  similarity  term  5,  since  the  image  pattern  is  not  correctly  deformed  and 
the  simple  model  of  deformation  (local  translation  and  constant  bv)  is  no 
longer  valuable.  The  measurements  are  then  projected  back  into  rotational 
motions  and  meridional  motions.  In  Fig.  7  we  compare  the  extracted  ro¬ 
tation  velocity  with  the  curve  of  the  theoretical  rotation.  Each  bin  is  0.5 
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degree  per  day  wide  and  1  degree  in  latitude  high.  The  relative  dispersion 
around  the  theoretical  rotation  is  approximatively  ±  5  %.  This  value  will  to 
be  compared  with  the  discrepancy  that  we  find  on  the  real  sequences  (Sec. 
4). 

4  EIT  Sequence  Analysis 

Three  successive  images  are  observed  on  the  1998-05-03,  repectively  at  times 
21:00:21  (ii),  21:12:09  ( 1 2)  and  21:25:35  ( I3 )  in  the  CME  Watch  mode,  at  the 
wavelength  19.5  nm  (see  Fig.  1).  A  second  sequence  of  two  images  (Jj.,  J2) 
observed  on  the  1999-04-17,  repectively  at  times  00:12:10  (Ji)  and  00:24:42 
(J2).  We  run  the  algorithm  on  rebinned  images  with  resolution  512  x  512 
pixels  to  reduce  the  computation  time.  We  choose  5  quality  bins  that  segment 
the  solar  disk  into  5  sets  of  equal  number  of  pixels. 

4.1  Differential  Rotation  Analysis 

We  analyse  the  sequence  (Jj,  J2)  in  the  same  conditions  of  section  3.  Figure 
8  shows  the  density  of  pixels  that  have  the  highest  estimation  quality  (levels 
4  and  5),  according  to  their  latitude  and  their  velocity  rotation.  We  also  plot 
the  theoretical  rotation  (plain  line).  The  estimated  velocities  are  variable 
(up  to  ±  2  degrees  per  day). 

We  note  that  for  the  quality  level  4,  the  rotation  velocities  of  the  two 
couples  are  ~  10  %  lower  than  the  expected  curve. 

Figure  9  (left),  the  subsampled  velocity  is  displayed.  After  examining 
the  histograms  (Fig.  8),  we  retrieve  the  cluster  of  pixels  that  have  a  fast 
rotation  velocity  (latitude  range  :  30  to  40  degrees).  A  zoom  on  a  particular 
region  belonging  to  this  cluster,  with  high  quality  estimation  (level  5)  and 
high  rotation  velocity  (above  15  degrees  per  day)  is  displayed  on  the  right, 
and  for  which  we  choose  to  overplot  a  greater  density  of  arrows. 

4.2  Brightness  Variation  Maps 

Figure  10  shows  the  BV  maps  of  the  two  successive  couples  of  images  (ii,  i2) 
and  (I2,  I3),  limited  to  the  ondisc  pixels.  In  the  latter,  an  obvious  change  in 
brightness  is  observed.  In  this  example,  the  thresholds  are  obvious  to  choose 
in  order  to  detect  brightenings  from  the  histogram  (Fig.  11).  The  same 
for  dimmings  or  negative  brightness  variations.  This  will  be  made  using  a 
set  of  selected  sequences  containing  typical  eruptive  events.  In  the  sequence 
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(images  h,h),  we  have  some  precursors  in  the  BV  maps  have  been  identified 
in  the  active  regions. 

Figure  11  (left),  the  plain  line  shows  peaks  (intensity  greater  than  1) 
of  outstanding  brightness  variations.  These  peaks  are  hardly  visible  in  the 
difference  image,  because  of  the  motion  effect,  which  motivated  the  use  of 
this  technique  to  process  long  sequence  of  EIT  images  (typically  one  day  of 
observations)  and  extract  the  interesting  but  faint  brightness  variation  areas 
from  the  bv  histogram.  In  the  right  figure,  the  extracted  region  is  outlined 
over  the  difference  image. 


5  Conclusion 

Real  and  semi-artificial  EIT  sequences  have  been  analysed  in  terms  of  mo¬ 
tions  and  brightness  variation  over  the  full  solar  disc.  This  estimation  is 
dense  (one  estimation  for  each  pixel),  but  a  quality  level  is  assigned  to  each 
pixel.  When  this  index  is  high,  the  estimation  is  reliable.  A  multiscale  im¬ 
plementation  refines  the  estimation  from  coarse  scale  to  fine  scales.  For  each 
pixel,  when  the  quality  estimated  at  one  scale  increases  at  the  finer  scale, 
then  the  estimation  is  updated.  Thus  the  algorithm  has  a  preference  for 
coarse  scale  estimations,  that  are  kept  through  the  multiscale  refinement  un¬ 
less  the  quality  of  the  estimation  is  locally  improved.  At  coarse  scales,  the 
number  of  pixels  used  in  the  estimation  is  higher  than  at  finer  scales:  the 
estimation  is  more  robust.  In  the  velocity  maps,  and  for  high  quality  levels  (4 
and  5),  we  observe  accumulations  of  pixels  into  clusters  that  have  a  common 
rotation  velocity.  Using  the  histogram  (Fig.  8),  we  have  identified  pixels 
from  an  active  region  (figure  9).  We  now  want  to  know  how  variable  is  the 
rotation  velocity  of  the  active  region  over  half  a  rotation.  This  will  be  pos¬ 
sible  after  the  processing  of  longer  sequences  (several  days) .  The  variations 
around  the  rotation  velocity  could  be  interpreted  as  oscillations  around  the 
theoretical  and  global  differential  rotation.  The  systematic  analysis  the  EIT 
archive  will  help  us  confirm  this  result.  Furthermore  a  full  disc  analysis  of 
our  algorithm  enables  an  operational  detection  of  events  in  sequences  of  two 
coronal  EUV  images,  using  the  BV  maps.  In  the  sequence  of  1998-05-03,  a 
strong  brightening  covering  a  large  area  has  been  undoubtedly  singled  out 
from  the  regular  solar  activity.  On  the  other  hand,  this  event  is  classified 
as  an  EIT  wave.  The  analysis  of  the  histogram  of  the  BV  map  shows  peaks 
deviating  from  regular  BV  histograms.  We  will  use  this  criterion  to  extract 
outstanding  intensity  variations  in  long  sequence  of  EIT  images.  A  future 
objective  of  this  is  to  associate  motions  and  brightness  variations  with  coro¬ 
nal  events  defined  by  observers.  We  will  also  investigate  the  correlations  in 
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time  between  the  outstanding  estimations  (fast  motions,  strong  brightness 
variations),  and  their  localizations,  and  the  in-situ  measurements  (LASCO 
CMEs  and  in-situ  ACE  data)  and  other  space  weather  observations  relevant 
to  event  prediction.  This  work  leads  to  solar  physics  research  as  well  as  real 
time  space  weather  services. 


A  PDE  formulation  of  the  optical  flow 


The  PDE  approach  considers  the  image  intensity  I(x(t),t)  as  a  physical 
quantity  that  is  conserved  over  time,  where  x  =  ( x ,  y)T  is  the  spatial  position 
in  the  image  plane.  The  PDE  form  of  the  BCA  is: 


<U 

dt 


=  0. 


The  differentiation  of  I(x(t),t)  in  Eq.  (22)  gives  the  OFCE: 


v/.ff+|  =  o. 


(22) 

(23) 


where  v  =  (dx/dt,dy/dt)r  is  the  velocity  vector  of  the  apparent  motion. 

As  the  real  signals  are  discretized,  the  temporal  and  spatial  derivatives 
must  be  approximated.  We  use  hereafter  the  finite  difference  scheme  to 
compute  the  temporal  derivative  dl/dt  using  the  difference  image  61  =  h~h 
and  the  time  distance  St  between  I\  and  / 2: 


dl  I(x,t  +  At)  —  I(x,t) 

—  =  lim  - ~r - - - 1 

dt  At— »o  At 


61  .. 
~  —  +  o(6t 
dt 


d2I 
dt 2 


)• 


The  discretized  OFCE  is: 


V/ -  u+  —  =  €(x,6t), 


(24) 


where  51  —  h(x)  —  h(x)  =  9a*  51.  The  order  of  the  discretization  error 
(different  from  the  model  error)  is  (see  Press  et  ah,  1992): 

a2/ 

a*,**)  ■ 

In  this  first  approach,  the  spatial  derivatives  can  be  computed  from  either 
image  A  and/or  /2.  In  the  case  of  image  registration,  the  order  of  error  due 
to  the  truncation  of  higher  order  terms  is: 


£(x,  5x)  ~  SxTHj  Sx, 
where  Hj  is  the  Hessian  matrix  of  I. 
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a.  The  gaussian  objects  numbered  from  0  to  7.  b.  Diffference  Image. 


Figure  5:  Left:  image  showing  the  8  gaussian  objects  G,  with  the  additive 
white  noise.  Right:  the  difference  image. 


Figure  6:  Left:  Quality  level  lines.  The  letter  ^  represents  the  i  -level  line  of 
quality.  At  ql  (lowest  level,  dotted-dashed  line),  some  noisy  structures  are 
still  selected.  The  level  q 2  is  a  dotted  line,  while  the  plain  lines  represent  the 
superior  levels  93,^4  and  q5.  Right:  Velocity  estimation  for  a  given  quality 
threshold  q. 
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a.  Rrigid  rotation.  b.  Synthetic  differential  rotation. 

Figure  7:  Calibration  procedure  using  semi- artificial  EIT  sequences  on  which 
a  synthetic  rotation  has  been  applied.  The  grey  level  represents  the  density 
of  solar  disc  pixels.  Vertical  axis:  synodic  rotation  velocity  (in  degree  per 
day).  Horizontal  axis:  solar  latitude  (in  degree). 


Figure  8:  Density  of  solar  disc  pixels  as  a  function  of  synodic  rotation  velocity 
(ordinate)  and  latitude  (abciss)  for  quality  levels  4  and  5. 
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Figure  9:  Left:  full  disc  subsampled  velocity  field  ( vx,vy ).  Right:  zoom  on 
an  active  region  with  a  rotation  velocity  greater  than  the  expected  value  of 
the  differential  rotation. 
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c.  BV(/2,/3)  d.  h-h 


Figure  10:  The  two  BV  maps  for  the  sequences  {1 1,12)  and  {h,h)  and  the 
corresponding  difference  images.  The  grey  levels  have  the  same  scale  in  each 
image. 
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Figure  11:  Left.  Plain  line:  histogram  of  BV  for  (/2 ,  ^3)  -  Dotted  line:  his¬ 
togram  of  the  difference  image  ( I3  —  h).  Dashed  line:  histogram  of  BV  for 
(/i,/2).  Right:  superimposed  on  I3,  the  thresholded  bv  map  of  (h,h)  using 
the  threshold  1  suggested  by  the  peaks  in  the  BV  histogram  (plain  line  in  the 
left  figure).  In  this  histogram,  the  brightness  variations  that  are  not  due  to 
the  motion  effect  appear  clearly,  while  they  are  not  observed  in  the  histogram 
of  the  difference  image  (dotted  line). 
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